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Abstract 

In this paper we address the issue of minimal time optimal control of fedbatch reactor in presence of 
complex non monotonic kinetics, that can be typically characterized by the combination of two Haldane 
models. The optimal synthesis may present several singular arcs. Global optimal trajectory results are 
provided on the basis of a numerical approach that considers an approximation method with smooth 
control inputs. 
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1 Introduction 

Fed-batch bioreactors represent an important class of bioprocesses, mainly in the food industry (e.g. yeast 
production or wine making) and in the pharmaceutical industry (like the production of the vaccine against 
the Hepatitis B) but also e.g. for biopolymer applications (PHB). It is also very much involved in the field of 
enzyme production which has been developed over the past decade due to the recombinant ADN technology 
and via the use of filamentous micro-organisms. One of the key issues in the operation of fed-batch reactors 
is to optimize the process operation over a limited time period. A intensive research activity has been 
devoted to optimal control of (fed-batch) bioreactors mainly in the seventies and in the eighties (see e.g. 
[SI IS] ) - In this paper we address the issue of minimal time optimal control of fedbatch reactor in presence 
of complex non monotonic kinetics, characterized here by the combination of two Haldane models, aimed to 
emphasize the presence of parallel metabolic pathways to transform the limiting substrate S into the biomass 
B. For those problems, it is well knwon that the optimal synthesis is bang-bang with the possibility of a 
singular arc [7J. For combinaisons of several non- monotonic growths, the multiplicity of singular arcs reveals 
a issue for determining which singular arc is globally optimal. In this work, global optimal trajectory results 
are provided on the basis of a numerical approach that considers an approximate problem whose optimal 
feedback laws are smooth. 

The paper is organized as follows. In Section [31 we present the model and the hypotheses. In Section 
131 we study the field of extremals, providing trajectories as candidate optimal solutions. Section U presents 
our approximation procedure and Section [5] shows how it can be applied to our problem, illustrated on two 
examples. We end by a conclusion in Section |B] 

*A. Rapaport is with the Equipe-projet INRA-INRIA 'MERE' (Modelisation Et Ressources en Eau) 
t Honorary Research Director FNRS, Belgium. 
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2 The model 



We consider a reaction scheme where n several bio-reactions consuming the same substrate S and producing 
a biomass B may occur simultaneously: 

S + B + Ri — >B + B i = l,---,n (1) 

Each reaction i requires an additional resource Ri, that we shall assume to be non-limiting, and is character- 
ized by a specific growth function Under the hypothesis of non-limitation by the auxiliary resources, 
one can assume that each function /iti(-) depends on the concentration of the substrate S only. 

Assumption 1. The functions Hi{.) are smooth, positive away from zero and null at 0. 

Moreover, we shall assume that each bio-reaction (that could be associated to different metabolic path- 
ways) transforms the substrate into biomass with the same yield y. The time evolution of the concentrations 
S and B in a perfectly mixed reactor, operated in fed-batch, is given by the following dynamical system: 

1 - O 

S = -~y>(5)B + £(S fn -S') (2) 
n n 

B = ^(S)B-fB (3) 

»=i 

V = Q (4) 

where V € (0, V max ] is the volume of the liquid phase in the tank, Si n > the input concentration of 
substrate and Q the input flow rate. 

In many applications (including wastewater treatment processes), a typical objective is to reach in min- 
imal time a target: 

{ (S, B, V) e R 3 + | S < Sref and V = V max } (5) 

where S% n > S re f > via the manipulated variable Q £ [0, Q max \- 
Let us define: 



/'• 



(S) = £>(S) (6) 



i=i 



x - f < 7 » 

the system dynamics ©-(HI) are equivalent to the one of a single bio-reaction with specific growth rate fi(-) 
and unitary yield factor: 

S = -KS)X + Q(S in ~S) (8) 

X = (i(S)X-^B (9) 

V = Q (10) 

One can easily check, from equations ©, that the quantity M = V(X + S — S in ) is constant along the 
trajectories, and consequently that it depends on the initial condition (Sq, X ,Vo) only. By considering the 
number M = V~o(X + So — Si n ), the dynamics that can be defined in the (S, V")-plane as follows: 

S = -v(S)(M /V-S + S in ) + Q(S m -S) (11) 
V = Q (12) 
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with the target: 

T= [0,5 re/ ] x {V max } (13) 

This minimal time problem has already been studied in [7] for cases where the function /i(-) has at most one 
maximum. An extension with impulse control inputs has been developed in [3]. With the help of a clock 
form and Green's Theorem, Moreno has proved that: 

1. the "bang-bang" strategy (i.e. Q = Q max until V — V max and then Q = 0) is optimal for any 
monotonic growth function /«(•), 

2. the "singular arc" strategy (that consists in reaching and remaining at S = S as long as possible) 
is optimal for growth functions increasing when S < S and decreasing when S > S (and under the 
condition Si„ > S > S re f). 

The clock form is used in Moreno's proof to show the global optimality of these strategies, a technique 
originated from the former results of Miele [5] (see also [J]). For cases where the growth functions have 
more than one local maximum on the interval (0, Si n ), this argument can still be used but only for the local 
optimality of singular arcs. Unfortunately, one cannot directly deduce the global optimality of the singular 
arc strategy. 

In the present work, we shall consider the following assumption. 

Assumption 2. The function /i(-) has a finite number of local maxima on the interval [0,5j n ]. Moreover, 
the set 

M = {S G [0, Si n ] I /•*(■) is locally maximal at S} (14) 

is such that: 

cardA4 > 1 (15) 

Sref < S- = minAl < S + = maxM < S in (16) 



This assumption can be typically fulfilled with n = 2 and /xi(-), A^O) of the Haldane type: 

= Ki + S + SP/U (17) 



that admits a single maximum at vA^Xi- For instance, the sum of the two Haldane functions 



Q.1 + S+10S 2 ' r v ; 5 + S' + 0.5 5' 2 
has two local maxima (see FigHJ. 



3 Study of the extremals 

We shall consider initial conditions on the domain T> — [S re f,Si n ) x (0,V rnax ] for our study (it sounds 
realistic that the initial substrate concentration is between the input and the desired ones). 

We first assume that the maximal flow rate Q max is large enough for ensuring the local controllability of 
the dynamics on the domain T>. 

Assumption 3. 

Qmax > max fl(S) ( — + Vmax) (18) 

se[S-,S + ] \b in -b J 
We show now that the target T can be replaced by a punctual one. 



3 



0.35 
0.30 
0.25 
0.20 
0.15 
0.10 
0.05 
0.00 




M-2 



M-i 



Figure 1: Graph of the sum of two Haldane functions. 



Proposition 1. From any initial condition in V, the optimal trajectory reaches the target T at point 



Proof. Consider the curve in the (S, F)-planc: 

C = {(a(V),V)\V e [0,V max ]} 
where cr(-) is solution of the differential equation 

MS) 




*>ref 



v Iv 



Then the domain £ delimited by this curve and the (S, V) axes: 

£ = {{S, V) | < S < a(V) and min{w > | a(v) >0}<V< V max } 



(19) 



(20) 



(21) 



is invariant whatever is the control. Assumptions 2 and 3 imply that the function er(-) is increasing and 
consequently S < S re f for any (S,V) in £. Note also that Assumption 2 imply that the function fi(-) is 
increasing on [0, S re f]. From the result of Moreno, one deduces that the curve C is an optimal trajectory. 

Consider an initial condition in T> \ T. It does not belong to £, but if an optimal trajectory reaches the 
target with S < Sref, it has to enter the domain £ i.e. it has to cross the curve C with S <C S re f , which 
contradicts the optimality of the curve C. □ 



Let us write the Hamiltonian of the optimization problem: 



H = A - \ s v{S)(M /V -S + S m ) + Q[X s 



Sir, 



V 



where Ao < 0, and the adjoint equations are: 

As = 
Ay = 

with 



As ( n'(S)X - n(S) + ® 



Ac 



-H(S)M + Q(S in - S) 



V 2 



X 



y i >->ir, 



(22) 

(23) 
(24) 

(25) 
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We define the switching function: 

= A s ^y^ + Ay (26) 
from which one deduces the optimality of "bang-bang" control inputs: 

if <b < 



max 



We focus now on the characterization of singular arcs 
Proposition 2. The singular arcs are trajectories 

S(t) = s 



if^>0 W 



V(t) = Vih)^*-^ + (eMS)(t-ti) _ i] 

Sin ~~ S V / 



M ° f„»(S)(t-tt) A (28) 



for t £ [ti,t2], with ti > t\ > 0, V{t\) < V ma x, V(ta) < V^ TO aa:, where SeM and the control input is given 
by the feedback equation: 

M 



Q°(S,V)=fx(S)l-2-^s +V ) ( 29 ) 

Proof. A singular arc strategy may occur when the switching function is identically equal to zero on a time 
interval of positive measure. One can easily compute: 

= A s ^y^'(S)X (30) 

Note from equation ((lip that from any initial condition in T>, one has: 

S(t) < S m , Vi > (31) 

and consequently one has: 

X(t) = ^ - S(t) + S m >0, Vt>0. (32) 

Note that it is not possible to reach the target from any initial condition in T> \ T, with a constant control 
input Q = and Q = Q ma x- Consequently, <f> has to be equal to zero at a certain time. Note also from 
equation (|23p that the sign of As is constant or As is identically equal to 0. In this latter case, Ay is constant 
and has to be non-zero from the Maximum Principle. Yet then <f> — Ay cannot be equal take to zero. So As 
is never equal to zero, and when tf> = one can conclude from H = 0, where H is the Hamiltonian defined 
in (|22|) . that As has to be negative. 

We deduce from equation (|30p that a necessary condition for an extremal to be singular is to have: 

fj,'(S) = (33) 

The zero of //(•) being isolated, we deduce that the singular arcs are of the form: 

S(t) = S with n'(S) = (34) 

Along with this condition, one can easily compute: 

4>= Xs{S ^^ fi"(S)XS (35) 

Assumption 3 implies that for tf> ^ 0, one has: 

sign(S") = sign(0) (36) 
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Consequently, in the neighborhood of S, one has: 

sign(^) = -sign(//'(S))sign(» 



(37) 



From the classification of fold points for time optimal control in the plane [T], one obtains that a point 
(S, V, As, Ay) such that 4> = is: 

- elliptic when fi"(S) > 0, and the optimal trajectory in its neighborhood is bang-bang, 

- hyperbolic when fi"(S) < 0, and the optimal trajectory in its neighborhood can have a singular arc. 

Finally, a necessary condition for an extremal to be singular is to have n'(S) = and (J-"(S) < which 
amounts to have S G Ai. 

Having S(t) = S on a time interval [^1,^2] implies S = 0, and from equation one deduces the 
expression of the control given in (|29p . Then from equation (|12l) the variable V is solution of the ordinary 
differential equation 

V = ti3)(^rg + V) , (38) 



whose explicit solution is given in (|28|) . □ 

Proposition 2 gives only a local optimality result. Away from the set M. x (0, V max ), we know that the 
optimal control is either or Q max and can switch, but we do not know a priori 

- towards which singular arc, defined by the value of S G M., it is optimal to go? 

- if is it optimal to quit a singular arc for reaching another singular one? 
To address the global optimality, we consider now a numerical approach. 



4 An approximation procedure 

Solving numerically minimal time problems with dynamics that are affine w.r.t. to the control input is 
usually intricate when one does know a priori the switching surfaces. The shooting function based on the 
integration of the Hamiltonian system is usually not smooth when the optimal control is discontinuous [12] . 
We consider here a smoothing method based on a idea originally proposed in [llj . We first define a new 
control: 

ui = j^--le[-l,l] (39) 



and denote: 



s 

V 



Then one can note that the dynamics (fTTjl - p"!?)) can be rewritten as follows 

i = F(£) + G 1 (0«i, £(0) = z e2? 

where 

- f 4S){^-s + s in ) 



F(0 = 

= 





s tn -s 

V 

1 



G(0 



(40) 
(41) 

(42) 
(43) 



(i 



According to Proposition [3J we shall consider the punctual target defined by: 



z f 



S, 



ref 



(44) 



The vector field Gi(-) is nowhere equal to the zero vector. Consequently there exists another vector field 
G 2 (-) such that: 

fect(Gi(0,G 2 (0) =» 2 , V^el? (45) 
One can also require G 2 (-) to be bounded: 

||G 2 (OII <r<+oo , V£eR 2 . (46) 

Now we consider the augmented dynamics, with an additional input w 2 : 

i = F(&) + Gi&Jm + eG 2 (Qu 2 , m = ^ (47) 

with u\ + u\ < 1 and e ^ 0. The Hamiltonian of this new problem is equal to: 

H e = PQ + + p*Gi (0« + ^G 2 (0v (48) 



The adjoint vector p being never equal to the zero vector (by the Maximum Principle), Condition (|45l) implies 
that: 

P 4 Gi(0 

Consequently the Hamiltonian iJ e is uniquely maximized by the smooth control inputs: 



(49) 



A/[p t Gi(0] 2 + [ep t G 2 (0] 2 
^G 2 (Q 

Vb i Gi(0] 2 + [ep t G 2 (0] 2 



(50) 
(51) 



We show now that the optimal trajectories for the extended dynamics converges toward an optimal trajec- 
tory of the original problem. This convergence has been recently studied in [TU], but the proof we propose 
here is different and is based on differential inclusions. 



Proposition 3. Let e n be a monotonic sequence of numbers converging to 0, and a sequence of optimal 
trajectories for e = e n with the same initial condition zq and target Zf. Then any £(•) limit of a sub-sequence, 
also denoted £„, in the following sense: 

£«.(•) — > £(•) uniformly and — > «;(•) weakly in L 1 (52) 

is an optimal trajectory for the original problem. Furthermore, there exists at least one such sub-sequence. 

Proof. Recall from Filippov's Lemma that the set of solutions of the equations (|4T[1 for measurable control 
inputs is exactly the set of absolutely continuous solutions of the differential inclusion: 

£ G * e (0 = F(0 + |J (Gf!(0«i + eG 2 (0«a) (53) 

Note that the set-valued maps Vf^ are monotonic w.r.t. e in the following sense: 

e<e' => *40c* e (0 (54) 

Let us consider an initial condition zq in T> and a monotonic sequence of positive numbers e„ converging to 
zero. As one has \&_ c (-) = 4 , e (-)) we can consider a decreasing sequence e n without any loss of generality. 
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Let us denote by T and T n the minimal times to reach Zf, respectively for dynamics (j^Tj) and (H71) for e = e n . 
Property (|54[) implies that the sequence T„ is non decreasing and bounded from above by T. Consequently 
T n converges to a limit, denoted T, such that T <T. 

Consider now a sequence £ n (") of optimal trajectories for the minimal time problem with e = e n . These 
trajectories can be prolonged up to time T (taking any admissible control input on the time interval [T n , T]), 
and are uniformly bounded on [0,T]. According to Dunford- Pettis Theorem, there exists a sub-sequence, 
also denoted £„(•) such that £«(•) converges weakly to v(-) on [0,T]. Let us then define: 

£;{t) = z Q + f v(s)ds ,t G [0,f] (55) 
Jo 

By weak convergence one has £ n (") — > £(•) an d £n(") - ^ £(■) weakly on [0,T]. One has also: 

L e *o(£„) + e™r-B a.e. * e [0, f ] (56) 
By compactness of trajectories of perturbed differential inclusions (see [5] or [T3]), one obtains 

Ce*o(C>.e. te[0,f] (57) 

Finally, one has 

Cn(r„) = Z, (58) 
and from the uniform convergence and continuity of trajectories £n('); one obtains 

ecn = (59) 

Consequently £(■) is an optimal trajectory for the original problem and necessarily T — T. □ 

Corollary 3. If the original problem admits a unique optimal trajectory £(•), then any sequence of optimal 
trajectories £«(•) for e„ a monotonic sequence of numbers converging to 0, converges uniformly to and 
£„(•) converges weakly to f(-). 

Remark. It is shown in [lOj that the optimal control does not necessarily converge point-wise in presence 
of singular arcs, and may exhibit a chattering phenomenon. Our approach here is slightly different, as we 
already know the locus of singular arcs and as our aim is to address the issues raised at the end of Section [3] 
This explains why we focus on the approximation of the optimal trajectories instead of the optimal controls. 



5 Numerical approximation 

In this section, we present a methodology that consists in computing the field of extremals for the regular- 
ized problem, and deducing if possible some properties of the optimal trajectories, namely answers to the 
questions raised at the end of Section [3] We illustrate this approach on two examples of growth function /u(-). 



We solve backward in time the Hamiltonian dynamics associated to the regularized problem: 

p = - P t d 5 F(0-p t d i G 1 (0<^p)-e nP t d 6 G 2 (Z)u* 2 (£, p) 

from the terminal state (£,p) = (zf,Pf) with different non zero vector pf. Without any loss of generality, 
one can choose vectors pj of norm equal to one i.e.: 



Pf 



cos(a) 
sin(a) 



(60) 
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Taking discrete values of a in [0,27r), we plot the projections of each solution in the (S, y)-plane. Two 
situations may happen : 

- either we fill the domain T> with a set of trajectories that do not intersect. Then each trajectory is optimal 
for the approximated problem, and consequently is close from an optimal one of the original problem (see 
Example 1 below). 

- either some trajectories intersect in D and only the part before the intersection (in backward time) can be 
optimal. Nevertheless, this partial information might be enough to answer the questions raised at the end 
of Section [3] (see Example 2 below) . 

For the vector field G%{-), we have simply chosen a constant one (but other choices are possible): 



G 2 (0 



(61) 



When G2O) is a constant, the adjoint equations are independent of e. Then one can show, similarly to the 
original problem, that p$ and pv are respectively negative and positive, for extremals with £(■) in the domain 
T>. Consequently, we take values of a in the interval (it, 3tt/2) only. 

In both examples below, values of parameters of the problem are given in Table Q] 

Table 1: Numerical simulation parameters 



Sref 


Vrnax 


Sin 


y 


M 


Qmax 


0.1 


50 


10 


5 


170 
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Example 1. We first test the method for a growth function /i(-) with only one maximum reached at S (see 
Fig GJ, for which the optimal solution is known [7]. 

1^ 




Figure 2: Graph of a function /i(-) with one maximum. 
It consists in reaching the singular arc S = S an staying on this arc until reaching the boundary of the 
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domain (i.e. V = V max ), using the following feedback 



if S > S 

Q S (S,V) if S = S and V < V max 
if S < S and V < V m ax 



(62) 



Extremals for the augmented dynamics with e = 0.01 are plotted on Fig|3j One can see that the domain 
T> is filled by extremals without intersection. Consequently each extremal is an optimal trajectory for the 
augmented dynamics, and one can check that they are close from the optimal trajectories for the original 
problem, given by the feedback (|62"j) . 




Figure 3: Extremals for the e-problem with one maximum. 



Example 2. We consider here a growth function /i(-) with two local maxima S\ and §2 (see Figure 2]). 

Extremals for the augmented dynamics with e = 0.01 are plotted on FigJ3] In this case, one can see that 
some extremals intersect, and so one cannot conclude on their optimality on the sub-domain B depicted in 
gray in Figure |5] 

Nevertheless, for the original problem (i.e. for e = 0), we know that away from the singular arcs S = S±, 
or S = S2, the optimal control is either or Q m ax (and can switch). So an optimal trajectory starting from 
S € (Si, §2) has to go toward S — Si or S = S2 (unless it touches the boundary V — V ma x) but we are not 
able to decide a priori 

- if it is optimal to reach V — V ma x before reaching a singular arc, 

- if not, towards which singular arc it is optimal to go, 

- if it is optimal to stay on a singular arc until reaching V = V max (it might be better to leave one 
singular arc to go to the other one). 

On Fig|3J one can also distinguish two areas A\ (resp. A2), in the complementary domain of the sub- 
domain B, depicted on FigJBl such that that extremals do not intersect and are close from the trajectories 
given by the feedback (|62]l with S = S% resp. 82- This observation is important because it allows us to 
conjecture that the optimal trajectories for the original problem reach one of the singular arcs and stay on 
it until reaching V = V max . 
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Conjecture. The optimal solution of the problem with a growth function presenting two local maxima at 
Si and 5*2 is given by the feedback Q§(-) where S £ {Si, §2}- 

Finally, the curve / on Figure [6] has been determined numerically as the set of points for which using 
feedback Q§(-) with S = Si or S = §2 give exactly the same time for reaching the target. Then, on the left 
part of the domain delimited by this curve, we conjecture that the control QgA-) is optimal, and Q§ 2 {-) in 
the right one. Furthermore, one can observe that when e get close from 0, the boundary of the sub-domain 
B get close from the curve I. 

6 Conclusion 

A novel numerical method for the investigation of minimal time problems in the plane, that may present 
several singular arcs, has been proposed. It is based on an approximation with no singular arc and for 
which extremals can be computed straightforwardly. The method has been applied on the optimal control of 
fed-batch processes with non-monotonic growth function. When the field of extremals of the approximated 
problem has no intersection, the optimal synthesis of the original problem can be deduced. Otherwise, we 
show that the method brings insights on optimal trajectories of the original problem on sub-domains of the 
state space. 
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